哈密尔顿蒙特卡洛的英文全称是Hamiltonian Monte Carlo，简称HMC，相比于MCMC（Markov Chain Monte Carlo，马尔可夫链蒙特卡洛）算法，HMC是一种更为快速的采样方法。

在哈密尔顿蒙特卡洛中，通常需要理解两个概念，一是哈密尔顿，二是蒙特卡洛。其中，这里的哈密尔顿来源于动力学，并对应着哈密尔顿动力学（即Hamiltonian Dynamics），常用于描述物体的运动过程。介绍哈密尔顿蒙特卡洛之前，我们先来简单地看一下MCMC算法和哈密尔顿动力学。

MCMC算法

MCMC算法是贝叶斯学习中的一种常用方法，能够用于描述物理系统，如下是一个与热能和温度 有关的系统： $$ p(\boldsymbol{x}) \propto \exp \left(-\frac{U(\boldsymbol{x})}{T}\right) $$ 其中，类似的分布形式被称为Gibbs canonical distribution（ "Gibbs正则分布" ）， $p(\boldsymbol{x})$ 表 示一个系统处于状态 $\boldsymbol{x}$ 的概率，并在数值上依赖于状态的能量 $U(\boldsymbol{x})$ (即热能) 和相应的温度 $T$ 另外， $p(\boldsymbol{x})$ 也能用来刻画物体关于位置和位置的运动过程，而在贝叶斯学习中，这个分布常用来 定义模型的参数，以多元正态分布为例， $$ p(\boldsymbol{x}) \propto \exp \left(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T \Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right) $$ 其中， $U(\boldsymbol{x})=\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T \Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu}) ， T=1$ ，当然，对于很多诸如基于能量的分布 (energy-based distribution) 都可以很轻松地写成Gibbs正则分布.

HMC原理

哈密尔顿动力学是利用给定物体在时刻 $t$ 的位置 (location) $\boldsymbol{x}$ 和速度 (velocity) $\boldsymbol{v}$ 来描述其 运动过程，由于位置 $\boldsymbol{x}$ 和速度 $\boldsymbol{v}$ 分别对应着势能和动能，将势能和动能分别描述为关于位置 $\boldsymbol{x}$ 和速度 $\boldsymbol{v}$ 的函数，则势能为 $U(\boldsymbol{x})$ ，动能为 $K(\boldsymbol{v})$. 系统的总能量 (势能+动能，常量，进一步求偏导可得到哈密尔顿动力学的偏微分方程) 被定义为 $H(\boldsymbol{x}, \boldsymbol{v})=U(\boldsymbol{x})+K(\boldsymbol{v})$ 在这里， $H(\boldsymbol{x}, \boldsymbol{v})$ 也被称为哈密尔顿函数，因此，位置 $\boldsymbol{x}$ 和速度 $\boldsymbol{v}$ 有独立的canonical distribution (正则分布)，即 $p(\boldsymbol{x}, \boldsymbol{v}) \propto \exp \left(\frac{-H(\boldsymbol{x}, \boldsymbol{v})}{T}\right)$ $=\exp \left(\frac{-U(\boldsymbol{x})}{T}\right) \exp \left(\frac{-K(\boldsymbol{v})}{T}\right)$ $\propto p(\boldsymbol{x}) p(\boldsymbol{v})$ 其中，我们发现，利用随机变量 $\boldsymbol{x}$ 和 $\boldsymbol{v}$ 的分布即可对联合分布进行采样。

以速度 $\boldsymbol{v}$ 为例，由于正则分布中的随机变量 $\boldsymbol{x}$ 和 $\boldsymbol{v}$ 相互独立，在哈密尔顿动力学中，若选择标准 正态分布为速度 $\boldsymbol{v}$ 的基本假设，则动能对应的函数为 $$ K(\boldsymbol{v})=-\log p(\boldsymbol{v}) \propto \frac{\boldsymbol{v}^T \boldsymbol{v}}{2} $$ 类似地，我们可以定义势能函数为 $$ U(\boldsymbol{x})=-\log p(\boldsymbol{x}) $$

HMC步骤

任意取一个 $q$ 的初始状态 重复如下计算:

• 1) 依据正态分布采样得到和
• 2）以和为初值，根据Hamiltonian方程计算L步，得到
• 3) 计算的接受拒绝概率: $$ \alpha=\min \left[1, \exp \left(-H\left(q_2, p_2\right)+H\left(q_1, p_1\right)\right)\right] $$
• 4) 从[0,1]均匀分布抽样得到一随机数 $u$ ，如果 $u \leq \alpha$ ，则接受采样，即赋值 $q_2 \rightarrow q_1$ ，回到第一 步。 很容易证明该算法满足如下细致平稳条件等式: $$ \begin{aligned} & \frac{1}{Z} \exp \left(-H\left(q_1, p_1\right)\right) \min \left[1, \exp \left(-H\left(q_2, p_2\right)+H\left(q_1, p_1\right)\right)\right] \\ =& \frac{1}{Z} \exp \left(-H\left(q_2, p_2\right)\right) \min \left[1, \exp \left(-H\left(q_1, p_1\right)+H\left(q_2, p_2\right)\right)\right] \end{aligned} $$ 之所以要进行接受拒绝操作是因为由于离散计算导致Hamiltonian系统Markov链转移概率发生变化。

HMC 实例

在 HMC中，通过构造出来的哈密尔顿动力学能量函数，能够很好地用到随机采样过程中的位置信 息和速度信息。现假设位置服从多元正态分布，即 $\boldsymbol{x} \sim \mathcal{N}(\boldsymbol{\mu}, \Sigma)$ ，其中，令均值为 $\boldsymbol{\mu}=\left[\begin{array}{l}0 \\ 0\end{array}\right]$ 协方差矩阵为 $\Sigma=\left[\begin{array}{cc}1 & 0.8 \\ 0.8 & 1\end{array}\right]$. 为了从 $p(\boldsymbol{x})$ 中进行采样，我们需要先确定势能 $U(\boldsymbol{x})$ (表征位置信息) 和 $\frac{\partial U(\boldsymbol{x})}{\partial x_i}$ (表征速度信 息）的表达式，前面我们已经将势能函数的正则形式构造成了 $$ U(\boldsymbol{x})=-\log p(\boldsymbol{x}) $$ 在这里，可以将势能函数写成 $$ U(\boldsymbol{x})=\frac{1}{2} \boldsymbol{x}^T \Sigma^{-1} \boldsymbol{x} $$ 对其求偏导数，有 $$ \frac{\partial U(\boldsymbol{x})}{\partial x_i}=x_i $$ 利用这些势能和势能的偏导数，就能够完成采样过程。

参考资料：

• https://zhuanlan.zhihu.com/p/32315762